// -*- C++ -*-

/*  *** *** *** *** *** *** *** *** *** *** *** *** *** ***  */
/*  JACOBI ROW CON ALLOCAZIONE STATICA della matrice nel     */
/*  modulo genera_A tramite la NEW:                          */
/*  static double tmpA[N][N];                                */
/*  *** *** *** *** *** *** *** *** *** *** *** *** *** ***  */
#define N 40
#define eps 0.01

#define MATRICE_IN   "/tmp/injacobi0"
#define MATRICE_OUT  "/tmp/outjacobi0"


generic main()
{
  stream double[N][N] A;
  stream double[N][N] risultato;

  genera_A      (                        output_stream A);
  jacobi        (input_stream A          output_stream risultato);
  raccogli      (input_stream risultato);
} 

 
genera_A (output_stream double A[N][N]) {
  proc_genera_A (output_stream A);
}


proc proc_genera_A (output_stream double A[N][N])
inc<"stdio.h", "iostream">
$c++{
  static double tmpA[N][N];

  FILE *fp;
  fp = fopen(MATRICE_IN,"r+");
  if (fp == NULL) 
    std::cerr << "Errore aprendo il file " << MATRICE_IN << std::endl;
  else {
    fread(tmpA,sizeof(double),N*N,fp);
    fclose(fp);
    assist_out(A, tmpA);
  }
}c++$


/* pragma parallelism degree 1 in jacobi; */

parmod jacobi (input_stream double A[N][N] output_stream double risultato[N][N]) {

  topology array [i:N] Pv;
  attribute double S[N][N] scatter S[*i0][] onto Pv[i0];  
  attribute bool diff replicated;
  stream double ris[N]; 

  init {
    VP i {
      init_false(out diff);
    }
  } 

  do input_section {
    guard1: on , , A {
      distribution A[*k0][*k1] scatter to S[k0][k1]; 
    }
  } while (true)
  
  virtual_processors {
    calcola_jacobi (in guard1 out ris) {
      VP i=0 {
	do {
	  perdi_tempo();
	} while (reduce (diff, ||) == true);
	assist_out (ris, S[i][]);
      }
      VP i=N-1 {
	do {
	  perdi_tempo();
	} while (reduce (diff, ||) == true);
	assist_out (ris, S[i][]);
      }
      VP i=1..N-2 {
	do {
	  calcola (in S[i][], S[i-1][], S[i+1][] out S[i][], diff);
	} while (reduce (diff, ||) == true);
	assist_out (ris, S[i][]);
      }
    }
  }
  
  output_section {
    collects ris from ALL Pv[i] {
      static double risult[N][N];
      const double *el;
      AST_FOR_EACH(el) {
	for(int j=0;j<N;++j)
	  risult[i][j]=el[j];
      }
      assist_out(risultato,risult);
    }<>;
  }
}


proc perdi_tempo() $c{}c$


raccogli (input_stream double risultato[N][N])
inc<"stdio.h", "iostream">
$c++{
  FILE *fD;
  fD = fopen(MATRICE_OUT,"w");

  fwrite(risultato,sizeof(double),N*N,fD);
  fclose(fD);

  std::cerr << "ENDING with SUCCESS j2.ast" << std::endl;
}c++$


proc calcola (in double el[N], double up[N], double dn[N] out double s[N], bool diff)
inc<"math.h">
$c++{
  bool t=false;

  s[0]=el[0];
  for (int j=1; j<N-1; j++) {
    s[j] = 0.20 * (el[j] + el[j-1] + el[j+1] + up[j] + dn[j]);
    t = t || (fabs(s[j] - el[j]) > eps);
  }
  s[N-1]=el[N-1];
  diff = t;
}c++$


proc init_false(out bool b) 
$c++{
  b=false;
}c++$


